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Abstract. This article gives an introduction for mathematicians interested in numerical computations in 
algebraic geometry and number theory to some recent progress in algorithmic number theory, emphasising 
the key role of approximate computations with modular curves and their Jacobians. These approximations 
are done in polynomial time in the dimension and the required number of significant digits. We explain 
the main ideas of how the approximations are done, illustrating them with examples, and we sketch some 
applications in number theory. 
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1. Introduction 

The purpose of this article is to give an introduction to the main results of the book [BCEJMj 
and their generalization in the PhD thesis [Brulj and in |Bru2j . as well as some applications, and 
most of all to explain the essential role played by approximate computations. The intended reader 
is a mathematician interested in numerical computations in algebraic geometry or number theory. 

The results concern fast algorithms in number theory and more precisely, fast computation of 
Fourier coefficients of modular forms. These coefficients, with Ramanujan's r-function as a typical 
example, have deep arithmetic significance and are important in various areas of mathematics, from 
number theory and algebraic geometry to combinatorics and lattices. 

The fastest previously known algorithms for computing these Fourier coefficients took expo- 
nential time, except in some special cases. The case of elliptic curves (Schoof's algorithm) was at 
the birth of elliptic curve cryptography around 1985. The results mentioned above give an algo- 
rithm for computing coefficients of modular forms in polynomial time. For example, Ramanujan's 
r(p) with p a prime number can be computed in time bounded by a fixed power of logp. 

Such fast computation of Fourier coefficients is itself based on the main result of the book: the 
computation, in polynomial time, of Galois representations over finite fields attached to modular 
forms by the Langlands program. 

The computation of the Galois representations uses their realisation, following Shimura and 
Deligne, in the torsion subgroup of Jacobian varieties of modular curves. The main challenge is 
then to perform the necessary computations in time polynomial in the dimension of these nonlin- 
ear algebraic varieties. Exact computations involving systems of polynomial equations in many 
variables take exponential time. This is avoided by numerical approximations with a precision 
that suffices to derive exact results from them. Bounds for the required precision - in other 
words, bounds for the height of the rational numbers that describe the Galois representation to 
be computed - are obtained from Arakelov theory. 

This article is organised as follows. Sections [2] and [3] are concerned with numerical methods 
used in the context of complex algebraic curves and their Jacobian varieties. Sections [5] and 
describe how to get exact results about torsion points on modular curves using these numerical 
methods. Section[3]focuses on the genus 1 curve An while Section[5]deals with the general modular 
curve Xi. As an application, Section [S] gives two examples of fast computation of coefficients of 
modular forms: Ramanujan's r-function, and the classical sums of squares problem. 
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2. Algorithms for curves and Jacobians 

Let X be a connected, smooth, projective algebraic curve over the field C of complex numbers. 
The set X(C) of complex points of X is a Riemann surface. Let g be the genus of X and let 
(wfc)i<fe< 9 be a basis for the space of holomorphic differentials on X. We fix a point b G -^(C) 
and we denote by Yb the set of homotopy classes of paths on X(C) starting at b. The universal 
cover fb :Yb — > AT(C) maps every path to its end point. The fundamental group tti(X(C), b) C Yb 
is the subset of (homotopy classes of) closed paths. It acts on Yb, with quotient X(C). We have 
an integration map fa : Yb — > C 9 defined by 

fa{l) = (/ Wl,..., / Ulg). 

The image of 7Ti(X(C), 6) by 0b is a lattice A in C 9 . It is called the lattice of periods. It is a free 
Z-module of rank 2g. The quotient C 9 /A is a complex torus. It is the set of complex points J(C) 
on the Jacobian variety J of X. The integration map 0b : Yb — > C 9 induces a map between the 
quotients X(C) — > J(C). This map is a morphism of varieties X — > J. We call this morphism fa 
also. For every positive integer k we denote §\ : X k — > J the morphism that maps (Pi, . . . , Pfc) 
onto fa (Pi) + • • • + 4>b{Pk)- Since the image in J does not depend on the ordering on the points 
Pj, we write X^ k ' for the fc-th symmetric power of X. We note that X^ is the quotient of 
X k by the action of the symmetric group. It is a nonsingular variety. We define the morphism 
0(*O . X {k) _^ j that maps |p i; _ ^ p k y Qnt0 ^(p^ + . . . + ^(Pfc). For fc = 5 the map b ff) is 
birational and surjective. It is not an isomorphism unless g < 1. Its fibers are projective linear 
spaces, mostly (but not all) points. A degree g effective divisor P = Pi + ■ ■ ■ + P g is said to be 

non-special if the map <f>^ is a local diffeomorphism at P. Otherwise we say that P is special. This 
definition does not depend on the chosen origin b. The set of special effective degree g divisors is 
the singular locus of <j)^ ■ All these maps </>j^ are called Abel-Jacobi maps. In particular 

4 9) ({Pi, . . . , P g }) = Yl ( f 3 mod A, 

where we can integrate § b Pj u>k along any path between b and Pj, provided we keep the same 
path for all k. We can apply the Abel-Jacobi map to any divisor on X. We set 4>b(J2j e jPj) = 

e j < t'b (Pj)- We note that for degree zero divisors, the image does not depend on the origin b. 
A divisor is said to be principal if it is the divisor of a non-zero meromorphic function on X . Two 
divisors are said to be linearly equivalent when their difference is principal. Any principal divisor 
has degree zero. A degree zero divisor is principal if and only its image by fa is zero. So the set 
J(C) = C 9 /A of complex points on the Jacobian is canonically identified with the group Pic°(X) 
of linear equivalence classes of degree zero divisors on X. 

We now list important algorithmic problems related to the Abel-Jacobi map. We illustrate 
them on the simple example of the projective curve X with equation 

(2.0.1) Y 2 Z - YZ 2 = a: 3 - X 2 Z. 

This curve has genus 1. We write x = X/Z and y = Y/Z. The unique (up to a multiplicative 
constant) holomorphic differential on X is 

dx dy 



2y-l x(3x-2) 

We choose the point b — [0 : 1 : 0] as origin for the integration map. For every computational 
problem we shall consider, we will also explain what can be proven when A' is a modular curve Xt 
and i (therefore g) tends to infinity. The definition of the modular curve Xi is given in Section [5j 
See also textbooks |Di-Shl IStej where Xe is often denoted Xi (I) . 
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2.1. Computing the lattice of periods. Wc first need a basis for the singular homology group 
Hi(X(C),Z). If X is the genus one curve given by equation (|2.0.ip . such a basis can be deduced 
from the study of the degree two map x : X — >• P 1 that sends (x,y) onto x and [0 : 1 : 0] to oo. 
This map is ramified at oo and the three roots of 4x 3 — Ax 2 + 1. We lift a simple loop around oo 
and one of these three roots. We then lift a simple loop around oo and another root. We thus 
obtain two elements in Hi(X(C), Z) that form a basis for it. 

Integrating a differential along a path is easy. We express the differential in terms of local 
coordinates. We then reduce to integrating converging power series. We integrate term by term. 
In case X is the curve given in equation (|2.U.1|) , we obtain a basis (fii,!^) for the lattice A of 
periods where 

fii = 6.346046521397767108443973084, 

fl 2 = -3.173023260698883554221986542+1.4588166169384952293308896131 
These calculation are made e.g. using the ■PARI] system. 

>al=0 ; a2=-l ; a3=-l ; a4=0 ; a6=0 ; 

>X= [al , a2 , a3 , a4 , a6] ; X=ellinit (X) ; 

>X. omega 

[6.346046521397767108443973084, 

-3 . 173023260698883554221986542 + 1 . 458816616938495229330889613*1] 

When dealing with general modular curves, an explicit basis for both the singular homology and 
the de Rham cohomology is provided by the theory of Manin symbols Man, Me reH ICrel IFrei fStej . 
Computing (good approximations of) periods is then achieved in time polynomial in the genus and 
the required accuracy |Cou2j . The practical side is described in Bosl, §6.3]. Textbooks [Cohj . 
[Crel Chapter 3] give even faster techniques for genus 1 curves, but we shall not need them. 

2.2. Computing with divisor classes. A degree zero divisor class can be represented by a 
point in the torus C 9 /A = J(C). It can also be represented by a divisor of the form 

(2.2.1) Pl + ... + Pg - gb 

in this class. This latter representation is not always unique. It is however unique for most classes 
because <pf is birational. The addition problem in this context is the following: given two degree 
g effective divisors P = P\ + • • • + P g and Q = Qi + ■ ■ ■ + Q g , one would like to compute a 
degree g effective divisor R = R± + • • • + R g such that the divisor class of R — gb is the sum of 
the divisor classes of P — gb and Q — gb. So we look for g complex points R±, . . . , R g such that 
Pi + ■ • • + P g + Qi + ■ ■ ■ + Qg — 2gb is linearly equivalent to i?i + • • • + R g — gb. This is achieved 
using the Brill-Noether algorithm [Bri-Noc, Vol . This algorithm uses a complete linear space C 
of forms or functions. This space should have dimension > 2g + 1. For example, assuming g > 4, 
we may take for C the space of all holomorphic quadratic differential forms. We compute once 
for all a basis for this space. Then the Brill-Noether algorithm alternates several steps of two 
different natures. Sometimes we are given a form (function) and we want to compute its divisor. 
Sometimes we are given an effective divisor D and we want to compute a basis for the subspace 
C(—D) consisting of forms (functions) vanishing at this divisor. 

The first problem (finding zeros of a given form) can be reduced, using a convenient coordinate 
system, to the following problem: given a power series f(z) = J2 k>0 fkZ k with radius of conver- 
gence > 1, find approximations of its zeros in the disk D(Q, 1/2) with center and radius 1/2. 
It is clear (see |Coull §5.4]) that, for the purpose of finding zeros, one can replace f(z) by its 
truncation J2o<k<K f^z k at a not too large order K. We then reduce to the classical problem of 
computing zeros of polynomials. A survey of this problem is given in [Coull §5.3]. 

The second problem (finding the subspace of functions vanishing at given points) boils down 
to finding the kernel of the matrix having entries the values of the functions in the chosen basis 
of C at the given points. 
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The only difficulty then is to control the conditioning of these two problems. This is done in 
two steps. We first prove [Coull §5.4] that the zeros of a holomorphic function on a closed disk 
are well conditioned unless this function is small everywhere on this disk. We then prove [Cou2, 
§12.7] that the form we consider cannot be small everywhere on any of the charts we consider, 
unless it has very small coordinates in the chosen basis of C. 

The resulting algorithm for computing in the group of divisor classes of modular curves is 
polynomial time in the genus and the required direct accuracy |Cou2l Theorem 12.9.1]. By direct 
accuracy we mean that the error is measured in the target space of the integration map, namely 
the torus C 9 /A. Saying that the direct accuracy is bounded from above by e means that the 
returned divisor R' = R[ + ■ ■ ■ + R' g is such that 

MR' -R) = $ ) {R')-$ ) {R) 

is bounded from above by e for the maxnorm in C 9 . This does not necessarily imply that the Rj 
are close to the Rj. Indeed, in case R = R\ + ■ ■ ■ + R g is special, there exists a non-trivial linear 

pencil of divisors R' such that 4^\r') — (j^ (R). Controlling the distance between R and R' will 
only be possible in some cases. 

In the special case when X is the curve given by equation (|2.0.ip the map <jy^ : X — » J is an 
isomorphism because the genus is 1. Computing with divisor classes is then very simple and the 
Brill-Noether algorithm takes a simple form. The space C consists of all degree 1 homogeneous 
forms, and a basis for it is made of the three projective coordinates X, Y and Z. Given P and 
Q, one considers the unique projective line Ai through P and Q. In case P = Q we take Ai to 
be the tangent to X at P. The line Ai meets X at three points: P, Q and a third point that 
we call S. We consider the unique projective line A2 through S and the origin b. The line A2 
meets X at three points: b, S and a third point that we call R. On can easily check that P + Q is 
linearly equivalent to b + R or cquivalently P — b + Q — b is equivalent to R — b. The coordinates 
of R can be computed using very simple formulae |Sil[ Chapter III]. We illustrate this using the 
PARI] system. We call P the point [0:0: 1]. We first compute Q such that Q — b is linearly 
equivalent to 2(P — b). We write Q — b = 2(P — b) using the = symbol for linear equivalence. We 
then compute R such that R — b = P — b + Q — b = 3(P — b). We then compute S such that 
S-b=Q-b+R-b= 5(P- b). 

>P=[0,0] ; 
>Q=elladd(X,P,P) 
[1, 1] 

>R=elladd(X,P,Q) 
[1, 0] 

>S=elladd(X,Q,R) 
[0] 

The answer for S means that S is just the origin b = [0 : 1 : 0]. So the divisor P — b has order 5 
in the Picard group Pic(AT), the group of divisors modulo linear equivalence. 

2.3. The direct Jacobi problem. Given a divisor on X we want to compute its image by <pb 
in the complex torus J(C) = C 9 /A. It suffices to explain what to do when the divisor consists of 
a single point P. For every 1 < k < g we then have to compute J, u>k- So we integrate uou along 
any path from b to P. We split the chosen path in several pieces according to the various charts in 
our atlas for the Riemann surface X(C). On every chart, the differentials u>k can be expressed in 

terms of the local coordinate. We then reduce to computing integrals of the form J 2 f(z)dz where 
f(z) is holomorphic on the unit disk. Such an integral can be computed term by term. When X 
is a modular curve, we have a convenient system of charts and a basis for C consisting of forms 
having small coefficients in their expansions at every chart. There is long standing tradition with 
stating and proving bounds for these coefficients. It culminates with the so-called Ramanujan 
conjecture. This conjecture was proved by Deligne as a consequence of |Dell| and his proof of the 
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analog of the Riemann hypothesis in the Weil conjectures in [Del2 . In case X is the elliptic curve 
given by equation (|2.0.ip we take for P the point [0:0:1] and find that 

<Pb ( p ) = / w = 2.538418608559106843377589234 mod A. 

Jb 

This integral is computed using the PARIJ system. 

> ellpointtoz(X, [0,0] ) 

2 . 538418608559106843377589234 

We notice that 

(D (F) = ^i modA 

So 5(P — b) is a principal divisor as already observed at the end of section [2731 

2.4. The inverse Jacobi problem. At this point we have two different ways of representing 
a degree zero class of equivalence of divisors. We can be given a divisor in this class like the one 
in equation (|2.2.ip . Such a divisor will be called a reduced divisor. We can also be given a vector 
in C 9 modulo the lattice of periods A. It is of course very easy to compute with such vectors. 
We also have seen in section 12.21 how to compute with reduced divisors. So both representations 
are convenient for computational purposes. We also have seen in section 12.31 how to pass from 
a reduced divisor to the corresponding point in the torus C 9 /A applying the Abel- Jacobi map. 
We now consider the inverse problem: given a point a mod A in the torus C 9 /A, find some 
P = Pi + ■ ■ ■ + P g such that the reduced divisor P — gb is mapped onto a mod A by <j)b- 

Using an iterative method We can try an iterative method like the secant's method. We 
illustrate the secant's method in case X is the curve given by equation (|2.0.1|) and 

(2.4.1) a = (Qi + n 2 )/ll = 0.2884566600635348685656+ 0.1326196924489541117573i. 

Starting from P = (50 - 50i, -223.147 + 547.739i) and P 1 = (20 - 20i, -54. 587 + 137.965i) we 
obtain an approximation up to 10 -26 after eighteen iterations. We use the |PARI| system and 
declare a function for the secant method. 

>secant (alpha, P0, PI, K)= 
{ 

local(f0,f l,xO,xl,x2,P2,P3) ; 
for(k=l,K, 

f0=ellpointtoz(X,P0) -alpha ;fl=ellpointtoz(X, PI) -alpha; 
x0=P0[l] ;xl=Pl[l] ; 
x2=xl-f l*(xl-x0)/(f 1-fO) ; 

P2=[x2,ellordinate(X,x2) [1]] ;P3=[x2,ellordinate(X,x2) [2]] ; 
if (abs (P2 [2] -P0 [2] ) > abs (P3 [2] -P0 [2] ) , P2=P3 , ) ; 
P0=P1;P1=P2; 

); 

return (P2) ; 
} 

The four parameters of this function are the target point in C/A, the two initial approximate 
values of P, and the number of iterations. We then type 

>alpha=(omegal+omega2)/ll ; 
x0=50-50*I;xl=20-20*I; 

P0= [x0,ellordinate(X,x0) [2]] ;Pl=[xl,ellordinate(X,xl) [2]] ; 
secant (alpha , PO , PI , 18) 
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Below are the results of iterations 14 to 18. We only give the values taken by the ^-coordinate. 

6.796891402429021881380876803 - 7 . 525836023544396684018482041i 

6.796539495414535904114103146 - 7 . 525907619429540863361002543i 

6.796539142100022043003057330 - 7 . 525908029913269174706910680i 

6.796539142094915910541452272 - 7 . 525908029899464322147329306i 

6.796539142094915911068237206 - 7 . 525908029899464321854796862i 



The continuation method Iterative methods only work if the starting approximation is close 
enough to the actual solution. Such an initial approximation can be provided by the solution of a 
different though close inverse problem. Coming back to our example, we will start from any point 
on X. Say Pq = (0, 0). We compute the image ao mod A of Po by the integration map. We then 
choose any P_i that is close enough to Po. 

>P0=[0,0] ; 

alpha0=ellpointtoz(X,P0) ; 
Pml=[0.1,ellordinate(X,0.1) [2]] ; 

We now move slowly from ao to a. We set a\ = ao + 0.1(a — ao) and we solve the inverse problem 
for ai using the secant's method with initial values P-\ and Po. 

>Pl=secant(alphaO+0. 1* (alpha-alphaO) ,Pml, P0,5) 

[0 . 218773824415936734050679268 - . 0122309960881052801981765895*1 , 
. 0388323642082357612959944279 - . 00390018046133107189481433241*1] 

We now set «2 = «o + 0.2 (a — ao) and we solve the inverse problem for a-i using the secant's 
method with initial values Po and Pi . 

>P2=secant (alpha0+0 . 2* (alpha-alphaO) , P0 , PI , 5) 

[0 . 410237833586311839505201998 - . 0205989424813431290064696558*1 , 
0.111775424533436210193603161 - 0.00838376796781394064004855129*1] 

We continue until we reach a 

>P3=secant (alpha0+0 . 3* (alpha-alphaO) , PI , P2 , 5) ; 

P9=secant(alpha0+0.9*(alpha-alpha0) ,P7, P8,5) ; 
P10=secant (alpha, P8, P9.10) 

[6 . 796539142094915911068237205 - 7 . 525908029899464321854796861*1 , 
-8 . 056577776742775028742861296 + 30 . 05694612451787404370259256*1] 

This continuation method is very likely to succeed provided the integration map has a nice local 
behaviour all along the path from ao to a. This is how practical computations have been realised 
in [Boslj for modular curves. It is however difficult to prove that this method works because the 
integration map (j)^ has a singular locus as soon as g > 1, and we do not know how to provably 
and efficiently find a path from ao mod A to a mod A that keeps away from the singular locus. 



3. Provably solving the inverse Jacobi problem 

We have presented in section 12.41 a heuristic algorithm for the inverse Jacobi problem. This 
algorithm is based on continuation. It seems difficult to prove it however because that would 
require a good control on the singular locus of the Jacobi map. In this section we present the 
algorithm introduced in Cou2] . This algorithm only requires a good control of the Jacobi map 
locally at a chosen divisor in X^- 9 K This is a much weaker condition and it is satisfied for modular 
curves. An important feature of this algorithm is the use of fast exponentiation rather than 
continuation. The principle of fast exponentiation is recalled in section 13.11 The algorithm for 
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the inverse Jacobi problem itself is given in section 13.21 Section 13.31 sketches the proof of this 
algorithm. Proving in this context means proving the existence of a Turing machine that returns 
a correct answer in a given time. One has to prove both the correctness of the result and a bound 
for the running time. This bound here will be polynomial in the genus of the curve and the 
required accuracy of the result. 

3.1. Fast exponentiation in groups. Assume we are given a group G. The group law in G will 
be denoted multiplicatively. We assume that G is computational. This means that we know how 
to represent elements in G, how to compare two given elements, how to invert a given element, 
and how to multiply two given elements. 

The exponentiation problem in G is the following: we are given an element g in G and an 
integer e > 2, and we want to compute g e as an element in G. A first possibility would be to set 
a\ = g and to compute au = asfc-i x g for 2 < k < e. This requires e — 1 multiplications in G. It 
is well known, however, that we can do much better. We write the expansion of e in base 2, 

e= 

0<k<K 

and we set bo = g and bk = b\_ x for 1 < k < K . We then notice that 

0<k<K 

So we can compute g e at the expense of a constant times log e operations in G. The algorithm 
above is called fast exponentiation and it admits many variants and improvements [Gor] . Its first 
known occurrence dates back to Pihgala's Chandah-sutra (before -200). See [DatSin, 1,13]. 

3.2. Solving the Jacobi inverse problem by linear algebra. Recall that we have two dif- 
ferent ways of representing an equivalence class of divisors of degree zero: reduced divisors or 
classes in the torus C 9 /A. We have seen that both models are computational. The Abel- Jacobi 
map 4>i '■ —> C 9 /A is computational also. We want to invert it (although we know it is not 
quite injective). More precisely we assume we are given some a in C 9 and we look for a degree g 
effective divisor on X such that <k (P) = <fib(P — gb) = a mod A. It seems difficult to prove the 
heuristic methods given in section ^. 4l for this purpose. So we present here a variant for which we 
can give a proof, at least when A is a modular curve Xg. We illustrate this method in the case 
where X is the curve given in equation (|2.0.ip . We still aim at the a given in equation (|2.4.1|) . 

We need a non-special effective divisor Pq of degree g. Since g = 1 we can take any point on 
X. For example Pq = (0,0). We note that the affine coordinate x is a local parameter at Pq. We 
choose a small real number e. The smaller e the better the precision of the final result. Here we 
choose e = 0.0001. We consider two points Pi and Pi that are very close to Pq. The first point 
Pi is obtained by adding e to the x-coordinate of Pq. The second point P^ is obtained by adding 
ei to the x-coordinate of Pq. 

P0=[0,0] ; 

P1=[0. 0001, ellordinate(X, 0.0001) [2]] ; 
P2=[0. 0001*1, ellordinate(X, 0.0001*1) [2]] ; 

We now compute the image a.\ mod A of Pi — Pq by the Abel- Jacobi map. We also compute the 
image ct2 mod A of P% — Pq. We note that a\ mod A is very close to 6 C/A. This is because 
Pq and Pi are close. We assume that a± is the smallest complex number in its class modulo A. 
We make the same assumption for ai. Then a\ and a<i are two small complex numbers, and they 
form an M-basis of C. This is because the integration map (f>^ is a local diffeomorphism at Pq (or 
equivalently Pq is a non-special divisor) and e has been chosen small enough. 

alphal=ellpointtoz(X,Pl)-ellpointtoz(X,PO) ; 
alpha2=ellpointtoz (X , P2) -ellpointtoz (X , P0) -omegal-omega2 ; 
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Recall that our target in the torus C/A is a mod A where a is the complex number given in 
equation (|2.4.ip . So we compute the two real coordinates of a in the basis {pt\, a.2). 

>M= [real (alphal) , real(alpha2) ; imag(alphal) , imag(alpha2)] ; 
coord=M~ (-1) * [real (alpha) , imag (alpha)] ~ 

[-2884.566581407009845250155464, -1326.196933330853847302268151]- 

We deduce that a is very close to a' — — 2884ai — 1326a2. And the class a' mod A is the image by 
4>b of — 2884(Pi — Po) — 1326(P2 — Po)- The linear equivalence class of the latter divisor is therefore 
a good approximation for our problem. There remains to compute a reduced divisor P — gb in this 
class using the methods presented in section |2"T21 Since the integers 2884 and 1326 are rather big, 
we use the fast exponentiation algorithm presented in section f3. 11 

>coord=truncate (coord) 
[-2884, -1326]- 

>Dl=ellsub(X,Pl,P0) ;D2=ellsub(X,P2,P0) ; 

P=elladd (X , ellpow (X , Dl , coord [1] ) , ellpow (X , D2 , coord [2] ) ) 

[6 . 798693122986621316758396123 - 7 . 528977879167267357619566769*1 , 

-8 . 059779911380488392224788509 + 30 . 07437308400090422713306570*1] 

We now check that the image of P — Po by <f>b is close to a 

>ellpointtoz(P) 

0.2884000018811813146007079855 + 0.1325999988977252987328424662*1 
>alpha 

0.2884566600635348685656351402 + 0.1326196924489541117573536012*1 

For a better approximation we should start with a smaller e. 

3.3. Matter of proof. The main concern when proving the algorithm in section I3~!?l is to prove 
that we can find an initial divisor Po that is non-special. In fact we must guarantee a quantified 
version of this non-speciality condition. The differential of at Pq should be non singular and 
its norm should not be too small. We can prove that such a condition holds true for modular curves 
|Cou2[ §12.6.7] because we have a very sharp description of these curves in the neighbourhood 
of the points called cusps. As a consequence we prove |Cou21 Theorem 12.10.5] that the inverse 
Jacobi problem for modular curves can be solved in deterministic polynomial time in the genus 
and the required direct accuracy. Recall that direct accuracy means that the error is measured in 
the target space C s /A. The main difference between the algorithm in this section and the one 
in section 12.41 is that we only need here to control the local behaviour of <j>^ at Po while the 
algorithm in section f2 .41 requires that the map (f){ be non-singular above the whole path from ao 
to a. 

In some cases it will be desirable to control the inverse error that is the error on the output 
divisor P in X^ 9 \ This will be possible when we can prove that 0(f) : X<*) -> J is a local 
diffcomorphism at P (that is P is non-special). We will also need a lower bound for the norm of 
the differential of 4>^ at P. Such a lower bound can be provided by arithmetic. 

4. Computing torsion points I 

In this and the next section we will assume that X is a modular curve and £ a prime number. We 
will be interested in ^-torsion points in the torus J(C) = C 9 /A. A point 

a = a mod A 

is an ^-torsion point if and only if a lies in iA. So the ^-torsion subgroup of J(C) is 7A/A and it 
has cardinality P 9 . This group is also denoted J[l}. 
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Some of these torsion points carry important arithmetic information. The values taken by 
algebraic functions at these points generate interesting number fields. We want to compute these 
fields. In this section we will focus on a special case. We will assume that X is the genus 1 curve 
given in equation (|2.0.ip and £ = 11. A more general situation will be studied in the next section^ 
We notice that the curve in equation (|2.0.ip is indeed the modular curve known as X%i. Since X 
has genus 1, the map 4>b : X — ¥ J is an isomorphism mapping b = [0 : 1 : 0] onto the origin. So the 
affine coordinate x and y induce algebraic functions io^ _1 and j/o^ 1 on J. There are ll 2 = 121 
points of 11 torsion in J and is one of them. We will be interested in the values taken by x o (p^ 1 
at the remaining 120 points of 11-torsion. On can check that x o (jj^ 1 takes the same value at two 
opposite points. So there only remain 60 values of interest. These are algebraic numbers and they 
form a single orbit under the action of the Galois group Gal(Q/Q). So it is natural to consider 
their annihilating polynomial 

(4.0.1) H{T)= [] (T-x^-Ha))). 

0^aSJ[ll]/±l 

This is an irreducible polynomial in Q[T]. Computing such polynomials is a cornerstone in the 
algorithmic of modular forms and Galois representations. 

4.1. An algebraic approach. The polynomial in equation (|4.0.ip is known as the 11-th division 
polynomial ipn of the genus one curve X. For every k > 1 one can define the fc-th division 
polynomial ipk (T) to be the annihilating polynomial of the ^-coordinates of all non-zero fc-torsion 
points on X. These polynomials can be computed using recursion formulae [Eng[ Section 3.6] [Sill 
Exercise 3.7] that follow from the simple algebraic form of the addition law on X. Using these 
recursion formulae we find 

H(T) = T 60 - 20T 59 + 112T 58 + 1855T 57 + • • • + 1321T 4 - 181T 3 + 22T 2 - 2T + 1/11. 

So we have an efficient algebraic method to compute H(T). We will explain in section [5] why it 
seems difficult to us to generalize this algebraic method to curves of higher genus. 

4.2. Using complex approximations. In this section we compute complex approximations of 
the coefficients of H(T). We also explain how one can deduce the exact value of these coefficients 
from a sharp enough complex approximation. We have seen in sections 12.41 and [3] how to invert 
the map 4>b- Given a point a in the torus C 9 /A we can compute a complex approximation of some 
reduced divisor P a — gb such that <t>b(P a — bg) = a. Since here the genus is one, P a consists of a 
single point on A, and it is uniquely defined. In case a = (f^i + f^)/!! we already found that the 
x-coordinate x(P a ) of P a is 

6.796539142094915911068237206 - 7.525908029899464321854796862i 

up to an error of 10~ 27 . We let a run over the 60 elements in ( J[ll] — {0})/ ± 1 and compute the 
60 corresponding values of x(P a ) with the same accuracy. We then compute their sum and find 
it is equal to 20 up to an error of 10~ 25 . This suggests that the coefficient of T 59 in H(T) is —20. 
In order to turn this heuristic into a proof, we need some information about the coefficients of 
H(T). We know that these coefficients are rational numbers. We need an upper bound on their 
height. The height of a rational number is the maximum of the absolute values of its numerator 
and denominator. We explain in the next section 14.31 how a good approximation and a good 
bound on the height suffice to characterise and compute a rational number. In case X is the curve 
given in equation (|2.0.ip an upper bound on the height of the coefficients of H(T) can be proved 
by elementary means. For example we know that the denominator of these coefficients is either 
1 or 11. In case A is a modular curve, similar bounds will be necessary. These bounds have 
been proved by the second author in collaboration with de Jong in |Ed-Jol| and [Ed- Jo2] . using 
Arakelov theory and arithmetic geometry together with a result of Merkl in [McrklJ on upper 
bounds for Green functions. 



10 



Jean- Marc Couveignes and Bas Edixhoven 



All the coefficients of H(T) are computed in the same way. They are symmetric functions of 
the x(P a ), so we can compute sharp approximations for them. We deduce their exact values using 
an a priori bound on their height. 



4.3. Recovering a rational number from a good approximation. In the previous sec- 
tion 14.21 we claimed that a rational number x = a/b can be recovered from a sharp enough 
complex approximation, provided we have an a priori bound on the height of x. We recall that the 
height of a rational number a/b, with a and b integers that are relatively prime, is max{|a|, \b\}. 
The rational number x — a/b is known if we know an upper bound h for its height and an approx- 
imation y of it (in R, say), with \x — y\ < l/(2h 2 ). Indeed, if x' — a! lb' also has height at most h, 
and x 1 7^ x, then 



a 


a' 




ab' - ba' 


b 


~ V 




W 



> 



We also note that there are good algorithms to deduce x from such a pair of an approximation y 
and a bound h, for example, by using continued fractions, as we will now explain. 

In practice we will use rational approximations y of x. Every rational number y can be written 
uniquely as 



[ao, 



fti, 



= a + 



1 



oi 



where n £ Z>o, ao G Z, a$ G Z>o for all i > 0, and a n > 1 if n > 0. To find these a*, one defines 
ao := ly\ and puts n = if y = ao; otherwise, one puts y± := 1/ (y — ao) and a\ = \_yi\ and n = 1 if 
yi = aj., and so on. The rational numbers [ao, ai, . . . , aj] with < i < n are called the convergents 
of the continued fraction of y. Then one has the following well-known result (see Theorem 184 
from [Ha-Wrj ). 

Proposition 4.3.1. Let y be in Q, a and b in Z with b ^ 0, and 

1 



y 



< 



2b 2 ' 



Then a/b is a convergent of the continued fraction ofy. 



5. Computing torsion points II 

In this section we describe how we compute the fields of definition of certain torsion points in 
Jacobians of modular curves. We recommend UDi-ShJ to those who are interested in an introduction 
to the theory of modular forms. 

Let SL2(Z) denote the group of 2 by 2 matrices with coefficients in Z and with determinant 
one. It acts on the complex upper half plane H via fractional linear transformations 

(5.D (: JV*=?^. 



c 



The standard fundamental domain F for SL2 (Z) acting on H (see Figure [1} consists of the z with 
\z\ > 1 and |3ff(z)| < 1/2. It is not bounded, hence not compact. Viewing H as the open northern 
hemisphere in P 1 (C) = C U {00}, with boundary the equator P 1 (R), we see that the closure F of 
F in P 1 (C) is the union of F and the point 00. 

For every prime number £ we let denote the subset of SL2 (Z) consisting of the ( " \ ) with 
c, a — 1 and d — 1 divisible by I. Then Ti is a subgroup of SL 2 (Z), of index £ 2 — 1. We assume 
that I > 5 from now on. Then the action of Ti on H is free. Each z in H has a neighbourhood U 
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Figure 1. Standard fundamental domain F for SL2(Z) acting on H 



such that all for 7 in Ti are disjoint. The quotient r^\H is therefore a Riemann surface that 
we denote by Yg, and the quotient map H — > Y/> is a covering map, that is, each point y in Yg has 
an open neighbourhood U such that the inverse image of U in HI is the disjoint union of copies of 
U , indexed by the inverse image of y. 

The Riemann surface Yg is not compact. A fundamental domain Fg in H for 1^ can be gotten 
as the union of the 7F, where 7 ranges over a set of representatives of Tg\ SL2(Z) /{l, — 1}. Such a 
set consists of [£ 2 — l)/2 elements and it can easily be found. We can compactify Yg to a compact 
Riemann surface Xg by adding I— 1 points, called cusps, the points of P 1 (R) that lie in the closure 
of Fg in P 1 (C). These points lie in fact in P 1 (Q) and can easily be written down. All this leads 
to an explicit topological and analytic description of Xg. It is covered by coordinate disks around 
the cusps. For example, the function 

(5.2) q-.M^C, z^e 2lxl \ 

restricted to the set of z with *<s(z) > l/£, induces a coordinate on a disk in Xg around the cusp 00. 
Indeed, the image under q of this region is the punctured disk of radius e~ 27T ^ e around 0, and the 
cusp 00 fills the puncture. The genus gg of Xg is equal to (£ — 5)(£ — 7)/24. For I = 11 the genus 
is 1, and indeed, X\\ is the elliptic curve X\\ given by equation (|2.0.1[) . 

It is of course a miracle that such an analytically defined Riemann surface as In is defined 
over Q, that is, can be described as a curve in a projective space given by a equations with 
coefficients in Q. But this is true for all £, and it is explained as follows, for £ > 13. The theory 
of modular forms gives that the C- vector spaces tt 1 (Xg) of holomorphic differentials on Xg have 
bases consisting of 1-forms ui whose pullback to H is of the form (J2 n>1 a n q n )-(dq)/q with all a n 
in Z. Quotients of such uj and us' in ^(Xg) then provide sufficiently many rational functions on 
Xg to embed it into a projective space, such that the image is given by homogeneous polynomial 
equations with coefficients in Q. 

We let Jg denote the Jacobian variety of Xg. It is also defined over Q, as well as its group law. 
This means that the group law is described by quotients of polynomials with coefficients in Q. 
Therefore, for all P and Q in Jg and for each a in Aut(C), the automorphism group of the field C, 
we have a(P + Q) = cr(P) + cr(Q). For each integer m > 1 the kernel Jg[m] of the multiplication 
by m map is finite (it consists of m 29e elements) and preserved by the action of Aut(C). This 
implies that all P in Jg[m) have coordinates in the algebraic closure Q of Q in C, that is, for each 
rational function / on Jg that is defined over Q and has no pole at P, the value /(P) of / at P is 
in Q. The analytic description above of Xg gives us an analytic description of Jg. 

We are interested in certain subgroups Vg of the ^-torsion subgroup Jg [£] of Jg that are invariant 
under the Galois group Gal(Q/Q) and consist of £ 2 elements. These Vg can be described explicitly 
and efficiently in terms of certain operators called Hecke operators on the first homology group 
of Xg. The whole point is to understand them algebraically, with their Gal(Q/Q)-action. 

The subgroup Vg defines a commutative Q-algebra Ag of dimension £ 2 as Q-vector space, the 
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coordinate ring of Vg over Q. This algebra Ag consists of the functions /: Vg — > Q with the 
property that for all a in Gal(Q/Q) and all P in Vg we have f{a{P)) = cr(/(P)). Addition and 
multiplication are pointwise. Each jg in Ag with the property that the fg(P) are all distinct is a 
generator, and At is then given as Q[T]/(Hf e ), with 

- n ( T -M p )) 

P&Vt 

A direct approach for computing An or Hf e algebraically, as in Section 14.11 in the case of the 
division polynomial "0111 is very unlikely to succeed in time polynomial in I, because in the case 
of Vg one has to work with the algebraic variety Jg, whose dimension grows quadratically with £. 
Writing down polynomial equations with coefficients in Q for Jg and Vg is probably still possible, 
in time polynomial in £. But computing a Q-basis of Ag from the equations in a standard way 
uses Groebner basis methods, which, as far as we know, take time exponential or even worse in 
the number of variables, that is, exponential or worse in I. 

For this reason we replace, in [BCEJMJ, exact computations by approximations. There are 
then two problems to be dealt with. The first is to show that jg can be chosen so that the 
logarithm of the height of the coefficients of Hf t , that is, the number of digits of their numerator 
and denominator, does not grow faster than a power of I. This problem is solved in Ed- Jolj . 
Merkl and [ Ed- Jo2j . using arithmetic algebraic geometry and analysis on Riemann surfaces. The 
second problem is to show that for the same choice of fg, the values fg(P) at all P in Vg can be 
approximated in C with a precision of n digits in time polynomial in n + 1. This is done in |Cou2j . 
The chapters [Bosl] and |Bos2) contain real computations using the method of Section 12.41 for 
prime numbers £ <T&. 

Let us now explain how we choose fg (up to some technicalities; the precise setup is given 
in |Ed3[ §8.2]) and say some words about the approximation of the fe(P). Standard functions on 
Jacobian varieties such as Jg are theta functions. But a problem is that these are usually given 
as power series in gg variables, and as gg grows this can make the number of terms that must be 
evaluated for a sufficiently good approximation grow exponentially in £. In other words, we know 
no method to approximate their values fast enough (of course, it is not excluded that such methods 
do exist). Our solution is to transfer the problem from Jg to Xf l , via the Abel-Jacobi map. We 
choose hg a suitable non-constant rational function on Xg, defined over Q, of small degree and with 
small coefficients. Then we take as origin a suitable divisor of degree gg on Xg, defined over Q. 
This divisor is carefully chosen in |Ed3| to have the following property: for each PinVg there is a 
unique effective divisor Qp = Q p.i + ■ ■ ■ + Qp, ge on Xg, such that its image under the Abel-Jacobi 
map is P. Then we define fg(P) — hg(Qp : i) + • • • + hg(Qp, gt ). Rather magically, the problem of 
power series in many variables has disappeared. The function hg is locally given by a power series 
in one variable. We evaluate it at each Qp t i separately. The Abel-Jacobi map (see Section [2]) is 
given by a sum of gg integrals of <7£-tuples of holomorphic 1-forms in one variable. The analytic 
description above of Xg and Jg should make it clear that the Abel-Jacobi map and the function 
hg can be well approximated with standard tools. That means that the only remaining problem is 
the inversion of the Abel-Jacobi map, that is, the approximation of the divisors Qp, but that was 
discussed and solved in Sections [2] and [3J The main result obtained in [BCEJM] is the following 
theorem. 

Theorem 5.3. There is a deterministic algorithm that on input a prime number £ > 11 computes 
the Q-algebra Ag in time polynomial in £. 

A probabilistic algorithm for computing Ag is also given in [BCEJM]. It relies on p-adic 
approximations rather than complex approximations. In |Cou3j it is explained how such p-adic 
approximations can be computed efficiently. From a theoretical point of view, a probabilistic 
algorithm is not quite as satisfactory as a deterministic one. From a practical point of view, it is 
just as good. In our case the probabilistic algorithm has a simpler proof than the deterministic 
one. And Peter Bruin [Brull IBru2] has been able to generalize it to a much wider class of Vg- 
likc modular spaces. Finding a similar generalization for the deterministic algorithm is an open 
problem at this time. 
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6. Applications and open questions 

The main motivation for all the work done in [BCEJM] is the application in number theory to 
the fast computation of coefficients of modular forms. Instead of attempting to present this in 
the most general case we give two examples: Ramanujan's r-function, and powers of Jacobi's 
^-function. 

Recall that q : H — > C is the function z M> e 2mz . The discriminant modular form A is the 
holomorphic function on H given by the converging infinite product 

(6.1) A = qH(l-q n ) 2 \ 

The holomorphic function A has a power series expansion in q, 

(6.2) A = J2r(n)q n , 

n>l 

whose coefficients, which are integers, define Ramanujan's r-function. It can be shown that for all 
( a j ) in SL2 (Z) , and for all z in H, we have 

(6.3) A^)-^-*,). 

Functions / : H — > C that are given by a power series J2 n >i a ™(/)<7™ with this symmetry under the 
action of SL/2(Z) on H with the exponent 12 replaced by an integer k are called cuspidal modular 
forms of weight k on SL2(Z). The complex vector spaces S^SI^Z), k) of cuspidal modular forms 
of weight k are finite dimensional. The dimension grows roughly as fc/12. More precisely, for 
k < 12 the space S^SI^Z), k) is zero, and S^SI^Z), 12) is one-dimensional, generated by A. The 
fact that each g in GL2(Q) with det(g) > acts on H and normalises SL2(Z) up to finite index 
leads to operators Tk. g on the S^SI^Z), k). These operators are named after Hecke. For each 
integer n > 1 there is an operator Tfc in ; for n prime it is the one induced by the matrix ( q 1 ) and 
for general n it is a bit more complicated. As the space S^SI^Z), 12) is one-dimensional, each 
T\2, n acts on it as multiplication by a scalar. This scalar turns out to be the coefficient r(n) of 
q n in the power series of A. Well known relations between the Hecke operators imply relations 
between the r(n) that are summarised in the identity of Dirichlct series, for s in C with real part 
3?(s) large enough: 

(6.4) Yl r (™)^ s = - t ^ s +p u p~ 2s r 1 - 

n>l p 

Here the product is over all prime numbers, and both sides converge for K(s) > 13/2. In fact, it 
is a famous theorem of Deligne ( [Dell] and [Del2]) that for all primes p one has 

(6.5) \ T (p)\ <2 P 11 '\ 

as conjectured by Ramanujan. 

The identity of Dirichlet series shows that for n and m with greatest common divisor 1 we have 
T(nm) — r(n)r(m), and that for p prime and r > 2 we have T(p r ) = t(p)t(jj t ~ 1 ) — p 11 r(p r_2 ). 
Therefore, the computation of r(n) is reduced to that of the r(p) for p dividing n. We can now 
state one of the main theorems of [BCEJM]. 

Theorem 6.6. There is a deterministic algorithm that on input an integer n > 1 together with 
its factorisation into prime factors computes r(n) in time polynomial in logn. 

Before this result, the fastest known algorithms to compute r(n) took time exponential in 
logn. For example, if one computes the product in (|6.1[) up to order n by multiplying the necessary 
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factors, then one spends time at least linear in n. To prove the theorem, it suffices to show that for 
p prime r(p) can be computed in time polynomial in log p. This will be done using Theorem 15.31 

The fact that modular forms have an enormous amount of symmetry as in (|6.3|) is certainly 
powerful, but it does not suffice at this point. What is needed is Galois symmetry, which is also 
behind Dcligne's famous result mentioned above. A lot could be said on this, but this is not an 
appropriate place for that. 

In a nutshell: modular forms give elements in de Rham cohomology of complex algebraic 
varieties defined over Q, and the singular homology with torsion coefficients Z/£Z of those complex 
varieties can be defined algebraically (Grothendieck, Artin, Deligne) and therefore has an action 
by Gal(Q/Q). 

For example, A gives rise, for every prime integer £ > 11, to a certain subgroup Vt of the 
^-torsion of the Jacobian Ji of Xg. This subgroup Vt has cardinality £ 2 . For £ ^ p the image of 
r(p) in Z/£Z is determined by the action of Gal(Q/Q) on this Vg. The addition map VexVf^Vi 
induces a Q-algebra morphism called co-addition from Ag to Ag <g> At, that is, from Q[T]/(Hf t ) 
to Q[T 1 ,T 2 }/(H fe (T 1 ),H fe (T 2 )). Computing r(p) modulo £ (for p ^ £) is then done by reducing 
Ai with its co-addition modulo p and computing on this reduction Ai iP a certain relation between 
the co-addition and the Frobenius map that sends a in Ag^ v to a p , just as in Schoof's algorithm for 
elliptic curves (see Section 1.2 of |Edll §1-2]). For more detail the interested reader is referred to 
Section 2.4 of |Ed2l §2.4] and the references therein. The point is that this advanced machinery 
can actually be used for computing r{p) mod £ in time polynomial in \ogp and £. 

In order to recover the actual value of r(p) as an integer, we compute r{p) modulo several 
small primes £. If the product of these small primes is bigger than 4p 5,5 then we deduce r(p) using 
inequality (|6.5|) and Chinese remainder theorem [Cohl 1.3.3]. 

We now come to our second example: the classical question in how many ways a positive 
integer n can be written as a sum of d > 1 squares of integers. Let us write rd(n) for this number, 
that is, r,}(n) = #{x 6 Z d : x\ + ■ ■ ■ + x 2 = n}. Then rd(n) is the coefficient of q n in the formal 
power series 9d, with 



(6.7) 9 d = J2rd(n)q n =J2l xl+ '" +x2d = E «* ? E ^ ) = e l in Z M- 

n>0 x£Z d WeZ / \x d e1 ) 

The formal power series 9\ defines a holomorphic complex function on the complex upper half 
plane 9 : H — > C by viewing q as the function q: z i— e 2mz . Poisson's summation formula then 
shows that for all z G H we have 

(6.8) 0(-l/4z) = (-2iz) 1 / 2 9(z), 



where the square root is continuous and positive for z in i-]R > o. This functional equation for 9, 
together with the obvious one 9(z + 1) = 9(z), imply that 9 is a modular form of weight 1/2, and 
therefore that 9d (interpreted as a function on H) is a modular form of weight d/2. 

This fact is the origin of many results concerning the numbers rd(n). The famous explicit 
formulas for the rd{n) for even d up to 10 due to Jacobi, Eisenstein and Liouville (see Mil and 
Chapter 20 of |Ha-Wr| ) owe their existence to it. In order to state these formulas, let J2d\m denote 
summation over the positive divisors d of m, with the convention that there are no such d if m is 
not an integer, and let x ■ Z — » C be the map that sends n to if n is even, to 1 if n is of the form 
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4rn + 1 and to —1 if n is of the form 4m — 1. Then we have: 
r a (n) =4^x(d), 

d\n 

r 4 (V)=8^d+16 E <2, 

2\d\n 2\d\(n/2) 

r 6 {n) = l&Y J x{^)d 2 -AY J X{d)d\ 

d\n d\n 

r s (n) = 16 ^ d 3 -32 E d 3 + 256 E d 3 , 

d\n d\(n/2) <*|(n/4) 

M^E*(^ 4 + f E*©^ E 

rf|n d|n d£l*[i], \d\ 2 =n 

In the last formula, Z[i] is the set of Gaussian integers a + bi in C with a and b in Z. 

Using these formulas, the numbers rd[n) for d in {2,4,6,8,10} can be computed in time 
polynomial in logn, if n is given with its factorisation in prime numbers. This is not the case for 
the formulas that were found a bit later by Glaisher for rd(n) for some even d > 12. We give the 
formula that he found for d = 12, as interpreted by Ramanujan: 

(6.9) n 2 (n) = 8E^ 5 -512 E + 16a "' where E a " 9 " = 9 II ^ ~ 9 2 ™) 12 - 

d|ra n>l m>l 

Computing a„ by multiplying out the factors 1 — q 2m up to order n takes time at least linear in 
n, hence exponential in logn. We know of no direct way to compute the a n in time polynomial 
in logn, even if n is given with its factorisation. However, J2 n>1 a n q n is a modular form, and 
therefore we can compute a n in time polynomial in logn, if n is given with its factorisation, via the 
computation of Galois representations. The same is true for the rd{n) for all even d. The explicit 
formulas for d < 10 correspond precisely to the cases where the Galois representations that occur 
are of dimension one, whereas for d > 12 genuine two-dimensional Galois representations always 
occur, as proved by Ila Varma in her master's thesis [Var] . 

We conclude that from an algorithmic perspective the classical problem of computing the rd(n) 
for even d and n given with its factorisation into primes is solved for all even d. The question as 
to the existence of formulas has a negative answer, but for computations this does not matter. 

Open questions Finally, we should point out that the algorithms in theorems 15.31 and 16.61 
despite their polynomial time complexity, are not so practical at present. However, Bosnian's 
computation of the Vg associated with A for £ in {13, 17, 19} enabled him to further study Lehmer's 
conjecture on the values of r(n) modulo n. See Lygeros and Rozier |Ly-Ro| for a more classical 
experimental approach. A challenge for the near future is to design and implement a practical 
variant of these algorithms. 

References 

[Bosl] J. Bosman, Computations with modular forms and Galois representations. In Computational As- 
pects of Modular Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton 
University Press, Princeton, NJ, 2011, 129-157. 

[Bos2] J. Bosman, Polynomials for projective representations of level one forms. In Computational As- 
pects of Modular Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton 
University Press, Princeton, NJ, 2011, 159-172. 

[BCEJM] B. Edixhoven and J.-M. Couveignes editors, with contributions by J. Bosman, J.-M. Couveignes, 
B. Edixhoven, R. de Jong and F. Merkl, Computational Aspects of Modular Forms and Galois 
Representations. Annals of Mathematics Studies 176, Princeton University Press, Princeton, NJ, 
2011. 



1G 



Jean- Marc Couveignes and Bas Edixhoven 



[Bri-Noe] A. Brill and M. Noether, Uber die algebraischen Functioncn und ihre Anwendung in der Ge- 
ometrie. Mathematische Annalen. 7 (1874), 269-310. 

[Brul] P. Bruin, Modular curves, Arakelov theory, algorithmic applications. PhD-thesis, Leiden, 2010. 
Available on-line at: http://hdl.handle.net/1887/15915 

[Bru2] Peter Bruin, Computing coefficients of modular forms. Publications mathematiques de Besancon. 
(2011), 19-36. Available on-line at: http://pmb.univ-fcomte.fr/2011.html 

[Coh] H. Cohen, A course in computational algebraic number theory. Graduate Texts in Mathematics 138. 
Springer, Berlin, 1993. 

[Coul] J.-M. Couveignes, Computing complex zeros of polynomials and power series. In Computational 
Aspects of Modular Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton 
University Press, Princeton, NJ, 2011, 95-128. 

[Cou2] J.-M. Couveignes, Approximating V/ over the complex numbers. In Computational Aspects of 
Modular Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton University 
Press, Princeton, NJ, 2011, 257-336. 

[Cou3] J.-M. Couveignes, Computing Vf modulo p. In Computational Aspects of Modular Forms and 
Galois Representations. Annals of Mathematics Studies 176, Princeton University Press, Princeton, 
NJ, 2011, 337-370. 

[Crc] J.E. Cremona, Algorithms for modular elliptic curves. Cambridge University Press, London, 1997. 

[DatSin] B. Datta and A.N. Singh, History of Hindu Mathematics. Motilal Banarsi Das, Lahore, 1935. 

[Dell] P. Deligne, Formes modulaires et representations Z-adiques. Seminaire Bourbaki 355 (1969). 

[Del2] P. Deligne, La conjecture de Weil. I. Inst. Hautes Etudes Sci. Publ. Math. 43 (1974), 273-307. 

[Di-Sh] F. Diamond and J. Shurman, A first course in modular forms. GTM 228, Springer, Berlin, 2005. 

[Edl] B. Edixhoven, Introduction, main results, context. In Computational Aspects of Modular Forms and 
Galois Representations. Annals of Mathematics Studies 176, Princeton University Press, Princeton, 
NJ, 2011, 1-27. 

[Ed2] B. Edixhoven, Modular curves, modular forms, lattices, Galois representations. In Computational 
Aspects of Modular Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton 
University Press, Princeton, NJ, 2011, 29-68. 

[Ed3] B. Edixhoven, Description of Xi(5l). In Computational Aspects of Modular Forms and Galois Rep- 
resentations. Annals of Mathematics Studies 176, Princeton University Press, Princeton, NJ, 2011, 
173-185. 

[Ed-Jol] B. Edixhoven and R. de Jong, Applying Arakelov theory. In Computational Aspects of Modular 
Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton University Press, 
Princeton, NJ, 2011, 187-201. 

[Ed-Jo2] B. Edixhoven and R. de Jong, Bounds for Arakelov invariants of modular curves. In Computa- 
tional Aspects of Modular Forms and Galois Representations. Annals of Mathematics Studies 176, 
Princeton University Press, Princeton, NJ, 2011, 217-256. 

[Eng] A. Enge, Elliptic curves and their applications to cryptography, an introduction. Kluwer Academic, 
New York, 1999. 

[Fre] G. Frey and M. Miiller, Arithmetic of modular curves and applications. In On Artin's conjecture for 
odd 2-dimensional representations. Lecture Notes in Math. 1585, Springer, Berlin, 1994. 

[Gor] D. M. Gordon, A Survey of Fast Exponentiation Methods. J. Algorithms 27(1) (1998), 129-146. 

[Ha-Wr] G.H. Hardy and E.M. Wright, An introduction to the theory of numbers. Fifth edition. Clarendon 
Press, New York, 1979. 

[Ly-Ro] N. Lygeros and O. Rozier, A new solution to the equation r(p) = (mod p). J. Integer Seq. 13 
(2010), no. 7, Article 10.7.4, 11 pp. 

[Man] Y. Manin, Parabolic points and zeta function of modular curves. Math. USSR Izvestija 6 (1972), 
no. 1, 19-64. 

[Merel] L. Merel, Universal Fourier expansions of modular forms. In On Artin's conjecture for odd 2- 
dimensional representations. Lecture Notes in Math. 1585, Springer, Berlin, 1994, 59-94. 



Approximate computations with modular curves 



17 



[Merkl] F. Merkl, An upper bound for Green functions on Riemann surfaces. In Computational Aspects of 
Modular Forms and Galois Representations. Annals of Mathematics Studies 176, Princeton University 
Press, Princeton, NJ, 2011, 203-215. 

[Mil] S.C. Milne, Infinite families of exact sums of squares formulas, Jacobi elliptic functions, continued 
fractions, and Schur functions. Ramanujan J. 6 (2002), no. 1, 7-149. 

[PARI] C. Batut, K. Belabas, D. Bernardi, H. Cohen, and M. Olivier, User's guide to PARI/GP (version 
2.3.1). http : / /pari . math . u-bordeaux . f r . 

[Sil] J. Silverman, The arithmetic of elliptic curves. Lecture Notes in Math. 106, Springer, Berlin, 1986. 

[Ste] W.A. Stein, Modular forms, a computational approach. With an appendix by Paul E. Gunnells. 
Graduate Studies in Mathematics, 79. American Mathematical Society, Providence, RI, 2007. 

[Var] I. Varma, Finding elementary formulas for theta functions associated to even sums of squares. Indag. 
Math. (N.S.) 22 (2011), 12-26. 

[Vol] E. Volcheck, Computing in the Jacobian of a plane algebraic curve. In Algorithmic Number Theory 
Conference. Lecture Notes in Computer Sciences, volume 877, Springer 1994, 221-233, 

Jean-Marc Couveignes, Univ. Bordeaux, 1MB, UMR 5251, F-33400 Talence, France. CNRS, 1MB, 
UMR 5251, F-33400 Talence, France. INRIA, F-33400 Talence, France 
E-mail: Jean-Marc.Couveignes@math.u-bordeauxl.fr 

Mathematisch Instituut, Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, Nederland 
E-mail: edix@ mat h . leidenuniv.nl 



